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\  A  commonly  occurring  problem  is  that  of  minimizing  least 
quares  expressions  subject  to  restrictions  on  the  solution. 


Dykstra  (1983)  has  given  a  simple  algorithm  for  solving  these  types 
of  problems  when  the  constraint  region  can  be  expressed  as  a  finite 
intersection  of  closed,  convex  cones.  Here  it  is  shown  that  this 
algorithm  must  still  work  correctly  even  when  each  cone  is  allowed 
to  be  arbitrarily  translated  (as  long  as  the  intersection  is  non¬ 
empty).  This  allows  the  algorithm  to  be  applied  to  a  much  larger 
collection  of  problems  than  previously  indicated.  C 
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An  Algorithm  for  Least  Squares  Projections 


onto  the  Intersection  of  Shifted  Convex  Cones 


Richard  L.  Dykstra 


1.  Introduction. 

The  problem  of  obtaining  least  squares  projections  subject  to 
various  constraints  is  a  frequently  occurring  problem  in  many  areas. 

For  example,  the  area  of  isotonic  regression  usually  is  concerned  with 
obtaining  least  squares  vectors  which  must  satisfy  certain  partial  order 
restrictions  (see  Barlow,  Bartholomew,  Bremmer  and  Brunk  (1972)  and 
Robertson  and  Wright  (l98l)).  Another  example  concerns  finding  the 
closest  (least  squares)  convex  (concave)  function  through  a  set  of 
points  in  the  plane  (see  Hildreth  (195^)  and  Wu  (1982)). 

Many  times,  the  constraint  region  can  be  written  as  a  finite  inter¬ 
section  of  simpler  constraint  regions.  This  raises  the  possibility  of 
using  iterative  schemes  based  upon  the  projections  onto  the  simpler 
regions  for  solving  the  overall  problem. 

John  Von  Neumann  (1950)  has  shown  that  if  the  constraint  region 
is  an  intersection  of  two  subspaces,  cyclic  iterative  projections  onto 
the  individual  subspaces  must  converge  to  the  desired  projection. 

Norbert  Weiner  (1955)  independently  proved  a  version  of  this  theorem 


in  a  slightly  different  setting. 
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Dykstra  and  Robertson  (1982)  have  developed  an  iterative  procedure 
for  finding  projections  of  rectangular  arrays  onto  the  class  of  arrays 
with  nondecreasing  rows  and  columns  based  only  upon  one-dimensional 
smoothings.  Later,  Dykstra  (1983)  extended  this  approach  to  the  general 
framework  of  projections  onto  the  intersection  of  closed  convex  cones. 

This  procedure  is  based  upon  finding  projections  onto  the  individual  cones, 
and  reduces  to  Von  Neumann’s  and  Weiner's  method  when  the  cones  are 
also  subspaces. 

It  is  the  purpose  of  this  paper  to  show  that  Dykstra' s  (1983) 
algorithm  can  be  extended  to  work  for  projections  onto  a  finite  non¬ 
empty  intersection  of  shifted  (translated)  closed,  convex  cones.  In 

particular,  this  means  that  least  squares  projections  under  general 

n 

linear  constraints  of  the  form  £  a.x.  <  (=)  b  can  be  handled  by 

1  1  1 

Dykstra' s  algorithm  even  when  b  4  0.  This  clearly  follows  by  writing 
n 

(x:  T  a.x.  <  b)  as 

1  1  1  “ 

n 

{x:  l  a.x.  <_  0}  -  (-b/a.  ,0,. . .  ,0) . 

1  1  1 

Constraints  of  the  form 

i  x.  <  zi,  i  =  l,...,n 


where  y  and  z  are  fixed  vectors  fall  into  the  translated  cone 
framework  by  writing 
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as 


( { x :  x  >  0}  -  (-y  ,-y  ,. 
1  ~  -L  «= 


,-y  ) )  n  ({x:  X.  <  0} 
n  i  — 


Other  examples  of  cone  constraints  given  in  Dykstra  (1983)  can  "be  general¬ 
ized  to  translated  cone  constraints. 


2.  Notation  and  Setting. 

We  denote  n-dimensional  real  coordinate  space  by  Rn ,  and  let 
g  and  w(w^  >0)  be  fixed  points  in  Rn.  The  inner  product  of  x 
and  y  (with  respect  to  w)  is  given  by 

n 

(x,y)  =  V  x.y.w  . 
i=l  111 


The  corresponding  inner  product  norm  of  x  is  defined  as 


(2.1) 


A  closed  (in  the  metric)  subset  K  of  Rn  is  a  closed  convex 
cone  if  x,y  K;  a,b  >_  0  implies  ax  +  by.  <s  K.  The  dual  cone  of  a 

closed,  convex  cone  K  is  defined  as 


(2.2)  K*  =  {y  Rn;  (y,x)  =  l  y.x.w  <0  Vx  <c  K>. 

1 

Of  course  K*  is  also  a  closed,  convex  cone  with  the  property  that 


K**  =  K. 


k 


A  commonly  occurring  problem  is  to  find  the  x  which  will 

(2.3)  Minimize  |jg  -  x|| 

x<£C 

where  C  is  a  closed,  convex  set  (x,y  C  C  implies  ax  +  (l-a)  <c  C 
Va  [0,1]).  A  vector  g*  5C  achieves  the  minimal  value  in  (2.3) 
iff 


(?.M  (g-g*,g*-f)  >  0  for  all  f5C. 

If  C  is  actually  a  closed,  convex  cone,  we  may  replace  (2.U)  by 

i )  ( g-g*,g *)  =  0,  and 

(2.5) 

ii)  (g-g*,f)  <  0  for  all  f  <s  C. 

(fee  Theorem  7.8  of  Barlow  et  al.  (1972).) 

Note  that  (2.5)  implies  that  if  C  is  a  closed,  convex  cone  and 
p*  solves  (2.3),  then  g-g*  e  C*. 


3.  The  Algorithm. 

We  wish  to  consider  problems  of  the  form 

(3.1)  Minimize  |j  g  -  x|| 

x  6  f  i  ) 
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where  K.  is  a  closed,  convex  cone  in  R  ,  b.  G  R  , 

1  1 

r 

K.-b.  =  {x-b.  ;  x  <cK.},  and  TifK.-b.)  ^  0.  We  assume  that  we  can 

ii  l  l  ^  l  l 

find  the  vector  in  K.-b.  which  will 

l  l 

(3.2)  Minimize  ||  f  -  x[| 

*€  K.-b. 

i  i 

for  any  f  and  any  i,  and  wish  to  use  this  ability  to  solve  (3.l). 

We  note  that  if  P(f  jc)  denotes  the  projection  of  f  Rn 
onto  the  closed,  convex  set  C,  then  for  a  closed,  convex  cone  K 

(3.3)  P(f  I  K.-b. )  =  P(f  +  b.  |K.)  -  b. . 

’ll  1  ’  1  1 

Thus 

(3.U)  f  -  P(f  |  Kj-bJ  =  (f  +  b.)  -  P(f  +  b.  |  K.)  <5  K*  for  all  f 
by  (2.5)  ii). 

We  shall  make  extensive  use  of  (3.**). 

Our  proposed  algorithm  is  identical  to  that  given  in  Dykstra  (19^3) 
except  that  we  allow  projections  onto  shifted,  closed ,  convex  cones. 

Our  scheme  can  be  succinctly  stated  with  the  aid  of  the  following 
notation: 

i)  For  any  positive  integer  n,  we  define  n(mod  r)  =  i  if 
n  =  kr  +  i  for  integers  k  and  i  where  1  <  i  <_  r. 
ii)  Initially,  set  n  =  0,  gQ  =  g,  and  1^  =  0  <s  Rn ,  i  =  l,...,r. 
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The  iterative  procedure  is  to 

1)  Set  Sn+1  =  P(gn  '  Tn(mod  r)  I  Kn(mod  r)  "  fcn(mod  r)}'  and 
then  update  Tn(mod  r)  by  resetting  it  equal  to 

°n+2  ^(mod  r )  ^ " 

2)  Replace  n  by  n+1  and  go  to  l). 

We  refer  the  reader  to  Dykstra  (1983)  for  further  elaboration  on  the 

algorithm  and  its  uses.  This  procedure  requires  only  the  ability  to  find 

projections  onto  the  (see  (3.3)).  These  individual  projections  are 

often  easy  to  program  and  quick  to  execute,  and  hence  can  be  combined 

to  solve  rather  difficult  optimization  problems. 

In  particular,  the  algorithm  applies  to  quadratic  programming 

n 

problem  with  a  finite  number  of  constraints  of  the  form  7  ax  <  b. 

1  1  1  " 

H.  Proof  of  the  Algorithm. 

To  simplify  the  proof  for  n  >  1,  we  will  write  g  in  (3.5) 

~  n 

aS  gk,i  When  n  =  (k-1)r  +  1  1  i  1  r.  (in  other  words,  .  is 

the  projection  onto  the  ith  shifted  cone  during  the  kth  cycle.) 

In  similar  fashion  we  depict 

*k,i  "  (gk,i-l~  if  2^^ 

^,1  "  ^^k-l,r  ~  Tk-1 ,1  ^  ’  lf  1  =  1* 

We  will  also  have  need  of  the  following  lemma,  a  proof  of  which  is  found 
in  Dykstra  (1983). 


7 


Lemma  U.l.  Suppose  a  sequence  of  nonnegative  real  numbers 
00 

r.  2 

is  such  that  )  a  <  <*>.  Then  there  exists  a  subsequence 
»  1  n 


(a  } 
n 


n=l 


U  > 

a  j=i 


such  that 


J 

Y  aa  -*-0  as  i 
t  m  n 
m=l  J 


We  now  establish  the  fundamental  result  of  the  paper. 


Theorem  U.l.  The  vectors  g  defined  in  (3.5)  converge  to  the 

n 

true  solution,  say  g*,  of  the  problem  defined  in  (3.1)  as  n  -*•  ®. 


r 

Proof.  Since  fl  (K.  -  b.)  4  0,  we  may  assume  WLOG  that 
r  1  1  X 

0  '  Pi  (K^  -  bj.  Then  b^t  and  the  true  solution  g*  exists 
uniquely.  Note  from  (U.l)  that 


(U  .2 ) 


,r ,  and 


gn-l,r  gn,l  In-l,l 


Thus ,  in  general  (I  .  =  0 ) ,  for  i  >_  2 

U  ,1 


IK. 1-1  -  8,11  2  ■  l|K.i-g*)  +  (I«-1.l-In.i,l|2 

*  2(8„,i *  V  Vi,i  -  - 2(e*  *  V  Vi.i  -  \,V 


Note  that  (g  .  +  b. ,  I  . )  =  0 
n,i  i  n,i 


(by  (2.5)). 


•8 


Moreover,  since  g  .  +  b,  £  K,  and  -I  ,  . .  €  K*  (by  (3.4)),  the 

n,i  1  1  n-l,i  i 

next  to  last  term  is  nonnegative. 

Similarly 


(4.4) 


ll g  ,  r  -  e*il2  i  llg„ ,  -  g*ll2  *  lli  ,  ,  - 1  ,ll2 

n-l  ,r  n,l  n-l,i  n,i 

-  2(g*  +  b  I  -  I  ). 

1  n-l ,  i  n  ,1 


Repeated  application  of  (4.3)  and  (4.4)  together  with  addition  and  the 
telescoping  property  of  the  last  term  yields 


n  r 


(4.5) 


ru  -  «*ir  >iiear-e*ir  +  r  i  \\\ !  t  -  \  j 

n’  k=i  i=i 


+2  7  ( g*  +  b . ,  I  .)  for  all  n . 

£  i  n,i 


Since  the  last  sum  is  nonnegative  (g*  +  b^.  €  ,  -1^  ^ K*),  we  know 


00  I* 


(4.6) 


Jx  ix111-,.-1.,.”  ' 


and  hence 


^n-l.i  '  Xn,ell  -  K.t-l  ’  ‘n.l11  Ui2)  and 


(4.7) 


I  ■‘'n-l,!  '  Tn  ,1  ^  =  ^6n-l,r  "  Sn,l^ 


■+  0  as  n  -*■  00  • 
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We  note  that  (U.5)  implies  that  g  and  g-g  =  I  ,+...+1 

n,r  n,r  n,l  n,J 

are  uniformly  hounded. 

Next  we  show,  that  there  exists  a  subsequence  {n  }  such  that 

J 


(U.8)  limsun  (I  .  +  I  0  +  . 

j  v1  v2 


.  +  I  ,  g  ,  -  f)  <  0  Vf  <c 

nrr  V1  “  i  1  1 


To  see  this,  note  that 


(I  +  •  •  •  +  I  ,  g  -  f ) 
n,l  n,r  n,l 


(U.9) 


Sn»1 


-f-(g  .+b.))( since (T  .  ,g  .  +  b. )  =  0 ) 


n  ,1  i 
r 


n,i  n,i 


=  j  (I  .  ,  g  .  -  g  . )  +  l  ( -I  . ,  f  +  b. ). 
■  i2  n  ,i  n  ,1  n  ,i  ^  n,i’  i 


The  last  sum  is  nonpositive  since  f  +  b.  G  K.  and  -I  ,  £k‘. 

i  l  n ,i  i 

For  the  first  part,  we  use  the  Cauchy-Schwarz  Inequality  to  say 


il2  <In,i>6n,l-8n,l)| 

i  J2  HT„,iH 

1  i=2  ' Jl  "I”'1  "  ‘  en,Jt1'  ’ 

=  (Jl  j2"Vi-  Vl.l  l|)(  j2llen,«-l  ■  “h.t11’ 


-  I 

m=l 


a  a 
m  n 


“V* 
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where  «n  =  JA,,-!  ‘  «n,JI  =  H \  ;£  '  Vl  ,£»  (s«"'-7))- 


Since 


a2  <  2r_2  l  ||l  -  I  || 2 

n  -  £“2  n,2  n-1 , 2 


(*+.6)  implies  that  £  a  <  °°.  Thus  lemma  ( 1+ .  1 )  can  be  employed  to 

n=l  n 

yield  (U.8).  Moreover,  since  the  ^  are  uniformly  bounded,  we  may 
assume  that  we  have  chosen  a  subsequence  such  that  (H.8)  holds  and 

gn  r  converges,  say  to  h.  Note  that  (**.?)  insures  that  g 

J’  n.  ,2 

J 

also  converges  to  h  for  every  2,  and  hence  that 
r 

h  <5  Pi  (K.  -  b.  )  (the  K.  are  closed).  In  addition,  since  the 
1  1  1  1 

g„  ™  -  g  =  I  _+...+  I  are  uniformly  bounded  and  converge 

Vr  v1  vr 

to  h-g,  we  may  use  ( 2.h )  to  argue  h  =  g*. 

Finally,  in  a  manner  similar  to  (U.5),  we  can  show 


(U.10) 


l«„  r  -  «*ll2  -  lie,  r  -  e*l|2  *  l  l  ||in  ,  -  i  ,  tl|2 

nj’  n’r  m=j+l  2=1  m  1,£ 

+  2  2  I  ((gm  2  +  V  “  (g*  +  V’  1  2  “  Tm  2}* 

m=n  +1  2=1  £  1  m~l’2  J"’5' 

J 


We  may  write  the  last  sum  as 


n  r 


( ^ • 11 )  2  l  l  (g  +  b  ,  I  ) 

m=nJ+l  2=1  m’£  *  m-1)£ 


n  r 


-2  * v  V’ 

m=nj+l  2=1  ’  ’ 
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(U.13)  +  2  1  (g*  +  b  I  ) 

o  x  n  9  A' 


(l4.ll*)  -  2  l  (g#  +  b  I  ). 

1=1  v 


Each  term  of  ( U .  11 )  Is  nonnegative  (g  .  +  b  <s  K„ ,  -T  <s  K*). 

in  y  %  Jt.  m-i  ^  *»>  ~ 

Each  term  of  ( U . 12 )  is  zero  by  (2.5).  Each  term  of  (U.13)  is  nonnegative 
for  the  same  reason  as  (U.ll). 

Finally,  we  can  write  (U.lU)  as 


-2[(g*,  l  I  )  +  l  (b  I  )] 
i=l  nj Z=1 


-2t(limg  2,  I  I  t)  *  l  (b  I  )]. 

j—  V  1=1  V  1=1  *  V 

-2[lim  [(g  ,  l  I  )  +  l  (b  ,1  )] 

j—  v1  *=1  V  1=1  V* 


-2  lim  [  l  (I  -,g  .  -  g  »)1 

n=2  v  v1  V* 


by  taking  f  to  be  Q  in  1*.9.  This  clearly  equals  zero  by  the 
00 

way  {n,}  was  chosen.  Then  noting  that  the  left  side  of  (1*.10) 

J  1 

goes  to  zero  as  J  it  easily  follows  that  g  -*■  g*  as  n  -*• 

n,r 

This  clearly  is  good  enough  by  (4.7). 
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